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PART  ONE 


THEORY  AND  EXAMPLES 

I.  STATEMENT  OF  THE  PROBLEM 

Figure  1 snows  the  problem  to  be  considered  and  defines  the  coordi- 
nates and  -tarameters  to  be  used.  The  perfectly  conducting  plate  covers 
the  entire  5=0  plane  except  for  the  aperture  which  Is  rectangular  In  shape 
with  side  lengths  L^Ax  and  L^Ay  In  the  x and  y d:'.rectlons  respectively. 

L and  L are  positive  Integers  and  L ^2.  The  aperture  Is  fed  by  a 
rectangular  waveguide.  The  excitation  of  the  waveguide  Is  a source  which 
produces  one  mode,  of  unit  amplitude,  which  travels  tov:ard  the  aperture. 

The  general  method  of  solution  [1]  Is  to  cover  the  aperture  with  a 
perfect  electric  conductor,  to  place  magnetic  current  sheets  and 
respectively  on  the  left-hand  and  right-hand  sides  of  this  conductor,  to 
obtain  an  Integral  equation  for  M by  equating  the  tangential  magnetic 
fields  on  both  sides  of  this  conductor,  and  to  solve  this  Integral  equa- 
tion using  the  method  of  moments.  The  testing  functions  are  the  same 
as  the  expansion  functions  for  M and  are  denoted  by  M^.  Each  ^ Is  a 
triangle  In  the  direction  of  current  flow  and  a pulse  in  the  direction 
perpendicular  to  current  flow. 

II.  SUMMARY  OF  BASIC  THEORY 

The  solution  for  M is  given  [IJ  by 

M = y V,M.  (1) 

<w  " 1/vw 

1 ' 

“V 

where  the  are  elements  of  a colun”  vector  V given  by 

^ = [y'^®  + . (2) 

Here, 

- - If  & • ■*" 

apert. 

[1]  R.  F.  Harrington  and  J.  R.  Mautz,  "A  Generalized  Network  Formulation 
for  Aperture  Problems,”  Scientific  Report  No.  8 on  Contract  F19628- 
73-C-0047  with  A.  F.  Caitbridge  Research  Laboratories,  Report  AFCRL- 
TR-75-0589,  November  1975. 
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A rectanj'iilar  waveguide  radiating  through  a rectangular 
aperture  into  half  space  bounded  by  a.i  electric  con- 
ductor. 


lis 

where  H ^(M.)  is  the  magnetic  field  at  radiated  by  the  magnetic 

M'] 

current  M.  in  front  of  the  perfect  electric  conductor.  Also, 


Ij 


I \ \ 


n 


in  n jn 


and 


where 


^^0  ^jO 


ij 


X e.  ds  . 


(4) 

(5) 

(6) 


apert. 


Here,  is  a unit  vector  in  the  z direction,  ^ is  the  jth  waveguide 

mode,  and  Y,  is  its  characteristic  admittance.  The  excited  mode  is  e„. 
j 

The  coefficients  of  the  reflected  waveguide  modes  are  given  by 


Wi 

(7) 

The  gain  associated  with  the  u component  of  the  magnetic  field  in  the 
half-space  z > 0 is  given  by 


Sttij  Re(V[Y^®]*^*) 


(8) 


where  n is  the  characteristic  Impedance  of  free  space  and 

P^“-2  ff  M^*ue  ^8  . (9) 

j JJ 

apert. 

In  (9) , points  from  the  distant  observation  point  to  the  aperture  and 
Ij^l  is  the  propagation  constant  k, 

[Y^®3  is  ^ times  the  admittance  matrix  dealt  with  in  [2].  has 

[2]  J.  R,  Mautz  and  R.  F.  Harrington,  "Electromagnetic  Transmission  Through 
a Rectangular  Aperture  in  a Perfectly  Conducting  Plane,"  Scientific 
Report  No.  10  on  Contract  F19628-73-C-0047  with  A.F.  Cambridge  Research 
Laboratories,  Report  AFCRL-TR-76-0056,  February  1976. 
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r> 


/'I 


been  calculated  in  [2]  for  the  special  case  = 0.  x^^  and  are 

defined  in  Fig.  1,  From  (9),  we  see  that  x^^  ^ 0,  y^^  0 merely  Introduces 

-jx,  (k  • u ) - jy.  (k  • u ) 

the  phase  factor  e which  la  subsequently  masked 

by  the  magnitude  operation  in  (8).  Now  that  [Y^®1  and  t are  disposed  of, 

m 

all  that  remains  to  carry  out  the  calculations  (1)  to  (9)  are  detailed 
formulas  for  A. , and  Y , . 

J 


III.  COEFFICIENTS  A. , AND  CHAKACTERISTIC  ADMITTANCES  Y. 
ij  J 

Expression  (6)  for  A^^  requt.:c3  a knowledge  of  the  expansion 


functions  1^  and  waveguide  modes  e^. 


The  set  1^  of  expansion  functions  is  split  into  a set  of  x 


^ ^ ^ 

directed  magnetic  currents  and  a set  ^ of  y directed  magnetic  currents 

defined  by 

^ V « fp=l>2, . . .L^-1 


V(q-l)(Vl)  “ ^ a=l  2 

X I q>=i,z » • • • 

M X/-.,  int  “ u T^(y->yJ  P (x-x.)  < 

-V  -I  1 P 1 1 

where  T^Cx)  and  T^(y)  are  triangle  functions  defined  by 


^ X - (p-l)Ax 


T^(x) 


Ax 

(p+l)Ax  - X 


Ax 

0 


(p-l)Ax  £ X pAx 
pAx  £ X £ (p+l)Ax 

|x  - pAxj  ^ Ax 


(10) 


(11) 


(12) 


ry  - (q-.l)Ay 
Ay 

Ty(y)  Jl.qti?.;y..--j^ 
q 1 Ay 

0 


(q-l)Ay  <.  y <.  qAy 

qAy  y <.  (q+i)Ay 
ly-qAy]  > Ay 


X V 

and  P (x)  and  P^(y)  are  pulse  functions  defined  by 
P q 


(13) 


4 


Pp(x)  =•( 


Pj(y) 


ri 

0 

1 

Lo 


(p-l)Ax  £ X < pAx 


all  other  x 
(q-l)Ay  £ y < qAy 


all  other  y 


(14) 


(15) 


The  set  e of  modes  for  the  rectangular  waveguide  is  split  Into 
TE 

a set  of  TE  modes  given  by  [3] 


TE 

^mfn(L  +1) 
m 


ab  e c 

m n r n mnx  . niry  m . mirx  mry, 

= X-  [u  :r  cos  sin  -u  — sin  cos  "r^] 

(mb)2  + <„a)2  t » b -y  a a b 

(16) 


m 


n ~ 0)  ly2y*«*X* 


m 


m + n 0 


n 

m = 0 

m = 1,2,, . . 


TM 


and  a set  e^  of  TM  modes  given  by 


TM  ^ - 

^nri-(n-l)L 

m 


ab 


(mb)  +(na) 
m =■  1,2,3, ...I 


r m raiTx  , mry  , n , mux  niryi 

7t[u  — cos  sin  —r^  + u — sin  cos  -r^J 

2W  a a b <"y  b a b 


(17) 


ra 


n “ 1,2,3, ...L 


n 


y Y TE  TM 

With  the  separation  of  into  ^ and  and  into  e^  and  , 

of  (6)  expands  to 


[3]  R.  F.  Harrington,  "Time-Harmonic  Electromagnetic  Fields,"  McGraw- 
Hill  Book  Company,  New  York,  1961,  Equations  (8-34),  (3-86),  and 
(3-89)  and  Section  4-3. 
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,xTE 


jx  „ TE  , 
M . • u X e , ds 

AV2  Mq 


apert< 


(18) 


yTE 


f • u 

J j ^ 


- ■<epds 

Mj 


apart. 


(19) 


.xTM 

'ij 


yTM 

ij 


ff 


TM  J 

_ X e.  ds 

MQ 


apart. 

‘ f 

apart. 


mV  TM  j 

M'  • u X e.  ds  . 

.«<L  <w<Z  V''^ 


(20) 


(21) 


Kara,  1^,  glvan  by  (10),  (11),  (16),  and  (17) 

raspactivaly. 

A tadlous  but  straightforward  avaluation  of  tha  intagrals  appearing 
in  (18)  to  (21)  leads  to 


/HH 

VAxAy  “Ij 


2u  .xTE 
A. 


ftTrAxAye 


ab 


ti  d niTT  V I 


/ sin 


mrAx  \ 2 
2a 


sin 


nirAy 

2b 


mirAx 

2a 


nirAy 

2b 


imt(x^  + pAx)  nii(y.  + (q-l/2)Ay) 
sin cos  : — r (22) 


i 


2ir  .yTE 
AxAy  ij 


f AirAxAye 


ab 


tn  . nil  X 
'“k.b^ 


, nirAy  \ 2 / , mirAx 
sin  \ f sin  -x— 


2b 


2a 


nirAy 

2b 


mirAx 

2a 


mii(x.,  + (p-l/2)Ax)  nir(y.  + qAy) 
cos sin  r (23) 
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I 2tt  .yTM  / 8ff AxAy  , mir .. 

yAxAy  ^Ij  "V  ab 


, nuAy  \ 2 / , irniAx  ' 

I irntAx 
V 2a 


rnrAy 

2b 


cos 


imr(x  + (p-l/2)Ax)  n-a(y.  + qAy) 

± sin  (25) 


where 


- -]/  (-)  + (f)  • 


(26) 


In  (22)  and  (24), 


1 = p + (q-l)(L  -1)  ■{ 


whereas  in  (23)  and  (25) , 


1 = p + (q-l)L 


In  (22)  and  (23), 


j “ in  + n (hjj+1) 


whereas  In  (24)  and  (25) , 


p « 1,2,...L  -1 


q “ 1,2 , . . .L 


r 


p = 1,2, ...L 


X 


q =•  l,2,...Ly-l 


m = 0,1,2, ...L 


m 


m+n  ^ 0 


n *”  0,1,2,***Xj 


n 


(27) 


(28) 


(29) 


j =■  m + (n-l)L 


m 


m •»  1,2, . . .L 


m 


n **  l,2,*«aXi 


n 


(30) 
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If  m or  n is  zero  in  (22)  or  (23) , the  resulting  (-■•^^--)  is  to  be  replaced 
by  unity. 


The  characteristic  admittni.ces  of  the  rectangular  waveguide  with 

relative  dielectric  constant  e and  relative  permeability  unity  are  classi- 

^TE  TM 

fied  as  either  TE  admittances  Y],  or  TM  admittances  Y given  by  [3] 


-jtlY 


TE 


j 


4)^  - 


f-^ 


. -W'r- 


j 


k >^e~  < k 
r J 


k > k. 
r j 


(31) 


r 


-jnifj  =•( 


-je 
**  r 

k,  - 

K 


k < k. 
r j 


(32) 


k /e”  > k.  . 
r j 


In  (31)  and  (32),  n is  the  characteristic  impedance  of  free  space,  k is 
the  free  space  wave  number,  and  k^  is  the  cutoff  wave  number  given  by  (26) . 
Strictly  speaking,  we  should  have  defined  separate  cutoff  wave  numbers, 
say  k^  and  k^  , for  TE  and  IM  modes  because  the  relationship  between  j, 

J J 

m,  and  n in  (26)  is  given  by  (29)  for  TE  modes  and  by  (30)  for  I'M  modes. 


IV.  POSSIBLE  DIFFICULTY  AT  CUTOFF  FREQUENCIES  OF  TM  MODES 

If  kj  =>  k*^,  then  Y™  of  (32)  tends  to  infinity.  Assume  that 
kj  is  very  close  to  k»^  for  j =»  j2>  j3>***jL  rewrite  [y'^®  + Y^®] 

of  (2)  as 

[Y^g  + yllBj  , jgg  4.  cj  (33) 

where 


8 


“ 1.2,3,,.. L 

(34) 

1 

C,.  + I Af  Y™  Af  + l'  A™  V™  Af 

Ij  Ij  ^ in  n jn  ^ In  n jn 

(35) 

i 

where  ^ denotes  the  sum  over  all  n except  n=j_,  j_,  j,,  ...  j . Also 
n 1 2 J L 

» 

1 

I 

.TE  .xTE  . , o /Y 

^In  “ '^in  * 1 “ 1,2,,..  (^^"l^Ly 

! 

[ 

1 

*S(V«L^,n  - C ’ i yW 

(36) 

! 

i 

IM 

and  Is  given  by  (36)  with  superscripts  TE  replaced  by  TM.  The 

^i^’  ^i™’  ^i™’  ^i™  defined  by  (18)  to  (21) 

Substituting  (33)  into  (2),  we  obtain 

1 

1 

1 

^ » [BB  + C]“^ 

(37) 

1 

Let 

^ - Ba  + ^2 

(38) 

and  choose  a such  that 

B - 0 

(39; 

i 

Next,  we  substitute  (38)  into  the  form 

- 

■=  [BS  + C]^ 

(40) 

1 

• 

of  (37)  and  use  (39)  to  obtain 

i 

i 

= [BB  + C]  Bo  + C^2 

(41) 

i 

i 

j 

! 

Expression  (41)  is  premultiplied  by  C"^  resulting  in 

! 

i 

V2  “ - [C"^BB  + U]  BO 

(42) 

I 


.9 


where  U Is  the  identity  matrix.  To  obtain  a,  premultiply  (42)  by  B and 
use  (39). 

a = [BB]**^[BC“^B  + U]'’^BC“^t^  ( 

The  solution  for  ^ is  obtained  by  substituting  (42)  and  (43)  into  (38). 

^ - C"^'B[BC“^B  + ( 


When  Y.  , i!.  = 1,2,... L,  tend  to  Infinity,  all  the  elements  of  B.  also 

h 

tend  to  infinity.  Hence  we  are  justified  in  substituting  the  approxi- 
mation 

[BC"^B  + U]“^  « [BC“^3“^  - [BC“^B]“^  + [BC“^B]"^  (45) 

where  the  superscripts  -2  and  -3  denote  respectively  the  square  and  cube 
of  the  Inverse  matrix,  into  (44)  with  the  result 

^ = C“^?^  - C"^B[BC"’^B]‘'^BC“^I^  + C‘"^B[BC“^B]"^SC’*^I^  - C"^B[BC”S]”^BC"^I^- 


Consider  the  case  in  which  the  waveguide  excitation  is  neither  the 

TM 

j-th  nor  the  j„th...  nor  the  j.th  TM  mode.  As  Y,  , A = 1,2,...L  approach 

1 2 L jj^ 

infinity,  (46)  approaches 


where 


Dij,  = aJ"  , &-l,2,3,...L 


ipu  TM 

Substitution  of  (47)  into  (7)  gives  for  all  j and  for  all  j 

except  j = j-,  jo  ...  j, . Here,  is  the  coefficient  of  the  jth 

TE  reflected  mode  and  is  the  coefficient  of  the  jth  TM  reflected 

J 

mode.  Expression  (47)  gives  T = 0,JI  = 1,2,...L.  To  investigate  the 

Jo 

TM  ^ 

manner  in  which  approaches  zero,  retain  one  more  term  in  (46). 
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+ C"^B[BC"^B]“^BC“^I^  (49) 

Premultiplying  (49)  by  B,  we  obtain 

B^  = [BC“^Br^BC“^t^  (50) 

which,  in  view  of  (7),  (34),  and  (48),  leads  to 


Y™  r™  - ([dc‘^d]"4c‘4^)  (5 

where  ( )^  denotes  the  £th  element  of  the  column  vector  inside  the 
parentheses.  From  (51)  we  conclude  that,  for  the  case  in  which  the 
waveguide  excitation  is  neither  the  j.th  nor  the  j.th  ...  nor  the 
jj^th  TM  mode,  as  Y ,£  = 1,2,...L  approach  infinity  the  magnetic  field 

of  the  i-  = 1,2,...L  TM  mode  approaches  a finite  limit  and  the  com- 

plex power  associated  with  this  mode  approaches  zero. 

Next,  consider  the  case  in  which  the  waveguide  excitation  is  the 
jpth,  1 1 p 1 L,  TM  mode.  From  (5)  and  (48),  is  2Y™  times  the  pth 

column  of  D.  The  column  vector  consisting  of  the  first  two  terms  in  (46) 

TM 

is  zero  because  it  is  2Y  times  the  pth  column  of  the  matrix 

- C"^D[DC”^]”^DC“^D 


all  of  whose  elements  are  precisely  zero. 
(46)  becomes 


TM 

Hence  for  large  Y.  , 1 

h 


*** 


^ = c“^B[BC“S]"^BC"4^  (52) 

which  reduces  to 

^ - 2(c"S[j)C"^D]"^)p  (53) 

where  ( ) denotes  the  pth  column  vector  of  the  matrix  inside  the 

parentheses.  Expression  (53)  gives  F]:  « 1 and  = 0,  i p.  By  retaln- 

Jp  Jjl 

Ing  the  last  term  in  (46)  we  obtain  the  more  accurate  representations 
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(54) 


,TM 


1 - ^ 

jm 


,TM 


where 


- 

yTM  ‘ 


i -p 


e,p  . 2 ([Sc-^i"^jp 


(55) 


(56) 


Here,  ( )j^p  denotes  the  Apth  element  of  the  matrix  Inside  the  parexi theses . 

From  (54)  and  (55),  v.  conclude  that,  for  the  case  In  which  the  wavegu.'de 

TM 

excitation  is  the  j th  TM  mode,  as  Y.  , S.  “ 1,2,...L  approach  infinity  the 
P 

aperture  electric  field  due  to  the  jj^th,  j2^h  ...  jj^th  TM  modes  approaches 
twice  the  electric  field  of  the  incident  mode  and  the  aperture  magnetic 
field  due  to  these  modes  approaches  a finite  limit.  The  complex  power 
associated  with  the  j^th  TM  waveguide  mode  approaches  23*^  and  that 
associated  with  the  j^th,  i p TM  waveguide  mode  approaches  zero.  The 
following  line  of  reasoning  shows  that  Real  (g  ) > 0.  From  power  con- 


PP 


sideratlons  [C  + C*]  is  positive^ definite . Hence  [C  ^ + [C  [Cc  '^D 

+ [0C~^D]*1,  and  t[0c  ^D]  ~ + [[DC  ^D]  ^1*]  are  in  turn  positive  definite 

and  finally  Real  (3  ) > 0. 

PP 

We  did  not  incorporate  any  of  the  logic  (33)  to  (56)  into  the  com- 
puter program  because  we  contend  that  the  simple  precaution  of  replacing 
R # 2 —6 

(-^)  - e of  (31)  and  (32)  by  10~°e  whenever  the  calculated  value  of 

Ic  IT  IT 


-i. 


e is  zero  is  adequ^ua.  Denoting  (-r^)  - e by  6 , we  justify 

t*  cc  r 3 


(%2  _ 

'•k'' 

this  contention  as  follows.  For  seven  decimal  digit  accuracy,  the  calcu- 
lated value  of  |6.(  must  be  exactly  zero  or  be  of  the  order  10  or 
J ,1  _5  ’• 

larger.  Hence,  we  assume  that  j6j|  is  roughly  10  or  larger.  In  the 

following  numerical  investigation,  what  appears  to  be  a reasonably  accurate 
solution  for  ^ was  obtained  for  |Sj|  in  the  neighborhood  of  10  We  set 


12 


L *L  = L = L “A,  Ax“Ay“  — ^ , a«b**-^>x,  "y,  “0 

K y m n » ^ ^ 


and 


= 1 such  that  **■  0 the  mode.  Successively,  we  let 


= -2.38  X 10"^,  -2.5  X lo"^,  - lo"^,  lo"^,  -lo"^^, 


-10 


-18 


The  first  value  -2.38  x lO  ^ was  the  spontaneously  calculated  value  of  6^. 
VJlien  the  waveguide  excitation  was  the  dominant  TE  mode,  the  elements  of 

■4* 

the  computed  solution  V for  the  first  four  values  of  6^  agreed  to  within  a 

fraction  of  a percent  of  the  magnitude  of  the  largest  element  of  In 
-12  ->■ 

going  to  6.  = - 10  , the  elements  of  V changed  by  only  a few  percent 

^ -18  -»■ 
whereas  for  6^  = - 10  , the  computed  solution  V was  absurdly  inaccurate. 

When  the  waveguide  excitation  was  the  TM  , mode,  the  first  four  values  of 

6 , gave  solutions  V differing  by  less  than  one  percent  whereas  the  com- 

^ *“12  ““18 
puted  solutions  V for  the  last  two  values  -10  and  -10  were  absurdly 

inaccurate. 


V.  SYMMETRY  OF  THE  ADMITTANCE  MATRIX 
The  admittance  matrix 

Y - [y''®  + Y^®] 

appearing  in  (2)  will  be  shown  to  be  symmetric. 

From  (A) , y'^®  Is  symmetric. 

llS  1 

Y ^ times  the  admittance  matrix  dealt  with  in  [2].  One  might 

argue  that  Y^®  is  symmetric  because  of  the  fact  that  the  set  of  testing 

functions  is  the  same  as  the  set  of  expansion  functions  and  because  of 

reciprocity.  However,  this  argument  is  not  strictly  correct  because  in 

[2]  the  integrations  over  source  and  field  regions  are  approximated  dif- 

hs 

ferently.  Nevertheless,  the  matrix  2Y  given  by  [2,  Eqs.  (23)  to  (26)] 
is  symmetric  because  one  can  show  that 


13 


XX  ^ XX 

ij  ji 


vyx  . 
ij  ji 


ll8 

The  symmetry  of  Y Is  not  due  to  reciprocity  but  due  to  the  fact  that 
of  two  rectangular  subareas,  the  first  as  seen  from  the  second  is  the 
same  as  the  second  as  seen  from  the  first. 


VI . CONTINUITY  OF  COMPLEX  POWER  FIOW 

It  will  be  shown  that  the  complex  power  flow  associated  with  the 
method  of  moments  solution  (2)  is  continuous  across  the  aperture. 

From  (2), 

- Iy”®  + Y*'®]^  (57) 

Premultiplying  the  complex  conjugate  of  (57)  by  we  obtain 

= v[y”^*  + Y^®*]^*  (58) 

which  is  the  same  as 

vt‘*  . vy”8*J*  . . (59) 

From  [1,  Eq.  (27)],  the  right-hand  side  of  (59)  is  the  complex  power 
radiated  into  half  space.  A development  similar  to  that  of  [1,  Eqs. 

(22)  to  (27)]  shows  that  the  left-hand  side  of  (59)  is  the  complex  power 
flowing  out  of  the  waveguide.  In  (59),  represents  the  interaction  of 

the  magnetic  current  with  the  Incident  magnetic  field  whereas  VY  * V repre- 
sents the  interaction  of  the  magnetic  current  with  the  magnetic  field  due 
to  the  magnetic  current. 

The  above  demonstration  of  continuity  of  complex  power  flow  is  based 
on  the  assumption  that  Y^®  is  given  by  (3)  and  that  y'^®  is  given  by  (4) . In 
reality,  the  integration  (3)  is  approximated  in  [2]  and  only  a finite  number 
of  teinns  of  the  infinite  sum  (4)  are  retained.  As  a result,  (59)  t-  still 
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true  but  the  right  and  left-hand  sides  of  (59)  are  merely  approximations 
to  the  complex  power  flow  on  both  sides  of  the  aperture.  This  is  tolerable 
because  the  obvious  way  to  obtain  more  accurate  complex  powers  would  be  to 
calculate  and  Y ® more  accurately.  Now,  these  more  accurate  y”®  and  Y^^ 
should  have  been  used  in  the  method  of  moments  solution  and  if  they  were  used 
then  the  right  and  left-hand  sides  of  (59)  would  be  more  accurate  approxima- 
tions to  the  complex  power  flow  on  both  sides  of  the  aperture. 


Since  the  left-hand  side  of  (59), 


henceforth  denoted  by  P 

wg 


P = op*  - vy”s*v* 

Wg 


(60) 


is  the  complex  power  flow  on  the  waveguide  side  of  the  aperture,  P should 

wg 


be  expressible  solely  in  terms  of  the  coefficients  of  the  reflected  wave- 
guide modes  and  the  characteristic  admittances  Y^ . Substituting  (4)  and  (5' 
into  (60) , we  obtain 

p = 2Y*  y A..V,  - y Y*  (y  a,  v.)  (y  a,  v*)  . (6d 

wg  0 ^ jO  j ^ n 'I  in  i'  jn  j" 


With  the  help  of  (7),  (61)  becomes 


r*ii2 

p = Y-(i  + rj  (1  - r.)  - y Y r 

wg  0 0'  ' 0 n'  n' 

which,  in  view  of  [1,  Eqs.  (22),  (45),  and  (48)],  is  not  surprising. 


(62) 


VII.  NUMBER  OF  WAVEGUIDE  MODES  REQUIRED 

The  nth  term  in  (4)  represents  the  contribution  of  the  nth  waveguide 

mode  to  y”®.  The  contributions  of  the  TE  and  TM  waveguide  modes  are 

mn  mn 

governed  by  the  (— - ■ — ) type  factors  appearing  in  (22)  to  (25).  Hence, 
these  contributions  begin  to  fade  away  gradually  as  soon  as 


mAx 

2a 


1 


2b 


15 


From  (44),  (34),  and  (7),  it  Is  evident  that  if  a solution  V gives 


r^-o. 


or 


1 + r 


j 


j 0 
0 , j - 


then  the  Jth  term  in  (4)  has  no  influence  on  ^ and  therefore  may  be 
neglected. 

If  the  aperture  is  centered  about  x “ a/2  or  y = b/2  or  both, 
then  it  will  be  shown  that  several  of  the  are  zero.  The  operator 
equation  [1,  Eq.  (4)] 

H^(M)  + ^(M)  - - hJ  (63) 

preserves  certain  symmetry  properties  of  the  transverse  component  of 
the  incident  magnetic  field.  For  Instance,  if  the  aperture  is  centered 
about  X = a/2  and  if  has  one  of  the  symmetry  properties 

1)  X component  even  about  x - a/2,  y component  odd  about  x = a/2 

2)  X component  odd  about  x •=  a/2,  y component  even  about  x = a/2 

then  the  magnetic  current  M has  the  same  symmetry  property.  Also,  if  the 
aperture  is  centered  about  y = b/2  and  if  h has  one  of  the  symmetry 
properties 

3)  X component  even  about  y ■ b/2,  y component  odd  about  y = b/2 

4)  X component  odd  about  y * b/2,  y component  even  about  y « b/2 


then  M has  the.  same  symmetry  property. 


The  above  symmetry  relations 

It  follows 


between  M and  will  be  verified  later  on  in  this  section. 

that  if  the  aperture  is  centered  about  x ■ a/2  and  if  IC  has  symmetry 

property  1)  then  T = 0 for  all  modes  which  have  symmetry  property  2)  whereas 
i 

if  has  symmetry  property  2)  then  F.  = 0 for  all  modes  which  have  symmetry 

property  1).  Likewise,  if  the  aperture  is  centered  about  y » b/2  then 

r.  = 0 for  all  modes  which  have  either  symmetry  property  4)  or  3)  depending 
J 1 i , 

on  whether  has  symmetry  property  3)  or  4) . Moreover,  because  H.  is  the 

transverse  magnetic  field  of  one  of  the  waveguide  modes  defined  by  [1,  Eq. 

4 

(45) ] , (16) , and  (17) , ^ has  either  symmetry  property  1)  or  2)  as  well  as 
either  symmetry  property  3)  or  4). 
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The  symmetry  relations  between  M and  will  now  be  verified.  As 
far  as  symmetry  properties  1),  2),  3),  and  4)  are  concerned,  has 

the  same  aymmetry  properties  as  M because  the  operator  merely  alters 
the  coefficients  of  the  expansion  of  ^ in  teirms  of  the  transverse  mag- 
netic fields  of  the  waveguide  modes  and  can  not  introduce  any  new  modes . 
From  [2,  Section  III], 


vAvere 


hNm)  . - 2J..F  - 27» 

-31' I Its' I 


F - 


apert 


[f 

|r-r' 

1.  II  I 

•F  JJ/  U-r'l 


ds 


apert. 


(64) 

(65) 


(66) 


P 


7^ 

-joi 


(67) 


where  r and  are  respectively  the  vectors  to  the  field  and  source 
points  in  the  aperture,  oi  is  the  angular  frequency,  e is  the  capacitivity 
of  free  space,  p is  the  permeability  of  free  space,  and  k “ wv^yF  is  the 
propagation  constant  in  free  space.  We  form  a table  of  symmetry  proper- 
ties of  p versus  those  of  M. 


SYMMETRY  PROPERTY  OF  M 

SYMMETRY  PROPERTY  OF  p 

1) 

odd  about  x ® a/2 

2) 

even  about  x = a/2 

3) 

even  about  y =•  b/2 

4) 

odd  about  y = b/2 

If  the  aperture  is  centered  about  x =>  a/2  then  iji  of  (66)  is  either  even 
or  odd  about  x >=  a/2  depending  on  whether  p is  even  or  odd  about  x =>  a/2. 
Similarly,  if  the  aperture  is  centered  about  y « b/2  then  (j>  is  either  even 
or  odd  about  y = b/2  depending  on  whether  p is  even  or  odd  about  y *=  b/2. 
Hence  <f>  has  the  same  symmetry  properties  as  p.  Tlierefore,  as  regards 
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symmetry  properties  1) , 2) , 3) , and  4) , V(|)  has  the  same  symmetry  proper- 
ties as  M.  Since  the  x and  y components  of  F^are  Integrals  of  the  same 
form  as  the  Integral  which  <|i  is,  F*u  has  the  same  even  or  odd  properties 

M5C 

about  X = a/2,  y *=  b/2  as  M*u  while  F»u  has  the  same  even  or  odd  proper- 

' ^ ^ My  r r 

ties  about  x = a/2,  y = b/2  as  M*u  . Thus  H^(M)  + H^(M)  has  the  same 

My  HA-t 

symmetry  properties  as  M.  From  this,  it  follows  that  M has  the  same  sym- 
metry properties  as 

VIII . SAMPLE  COMPUTATIONS 

A computer  program  using  the  formulas  derived  in  the  preceding 
sections  has  been  written.  It  is  described  and  listed  in  Part  II  of  this 
report.  In  this  section  we  give  some  examples  of  the  computations  that 
can  be  made  using  the  general  program. 

Figure  2 shows  the  equivalent  magnetic  current  for  a rectangular 
waveguide  of  dimensions  X by  X/2  radiating  into  half  space  through  a narrow 
centered  rectangular  slot  of  dimensions  X by  X/10.  Figure  2(a)  shows  the 
x-component  of  equivalent  magnetic  current,  which  is  also  equal  to  the  y- 
component  of  tangential  field  in  the  slot.  No  y-component  of  magnetic 
current  was  obtained  because  only  one  pulse  in  y was  used.  M is  normalized 
with  respect  to 


(68) 


where  the  Integral  is  over  the  waveguide  cross  section.  In  other  words, 
the  normalization  factor  is  the  root-mean-square  value  of  the  incident  E 
field.  The  phase  of  M is  with  respect  to  that  of  the  Incident  electric 
field  at  the  aperture.  All  computations  are  for  dominant  mode  ex- 
citation. Figure  2(b)  shows  the  radiation  gain  patterns  in  the  two  planes 

X “ 0 and  y ■ 0.  The  notation  denotes  the  gain  pattern  due  to  Hg  in 

the  y - 0 plane.  The  notation  G denotes  the  gain  pattern  due  to  H in 

XX  X 

the  X - 0 plane.  The  horizontal  axis  in  Figure  2 is  the  z axis. 
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(b)  0.8A 

Fig.  4.  Radiation  gain  patterns  for  an  open-ended 
rectangular  waveguide  of  dimensions  a/b  » 
2.25  radiating  into  half  space. 
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CROWI-EV.  LEVIS 
COMPUTER  PROGRAM 


The  equivalent  aperture  admittance  seen  by  the  dominant  mode 
for  an  open-ended  square  waveguide  of  width  a radiating  into 
half  space.  Our  computed  results  are  compared  to  those  cal- 
culated and  measured  by  Cohen,  Crowley,  and  Levis  [4], 


Figure  3 shows  the  equivalent  magnetic  current  for  an  open-ended 
rectangular  waveguide  of  dimensions  a/b  “2.25  radiating  into  half  space, 
for  two  different  wavelengtha.  Figure  3(a)  is  for  a ■ 0.6X,  and  Figure  3(b) 
is  for  a “ 0.8X.  Although  the  aperture  is  now  relatively  wide,  the  y- 
component  of  magnetic  current  is  still  small,  and  thus  is  not  shown.  The 
top  Figure  3(a)  shows  magnitude  and  phase  of  the  x-component  of  M at  y ■=  b/2, 
and  the  bottom  Figure  3(a)  shows  them  for  the  x-component  of  M at  x = a/2. 
The  top  Figure  3(b)  shows  the  magnitude  and  phase  of  the  x-com:  ant  of  M at 

y “ 3b/8,  and  the  bottom  Figure  3(b)  shows  them  for  the  x-component  of  ^ at 
X “ a/2.  Note  that  M is  zero  at  x “ 0 and  x ■ a,  and  it  is  large  near  both 
y “ 0 and  y - b.  This  is  to  be  expected  from  theory.  Again  the  magnetic 
current  is  normalized  according  to  (68) . Figure  4 shows  the  principal  plane 
radiation  patterns  for  the  same  waveguide-fed  apertures,  4(a)  for  a “ 0.6X 
and  4(b)  for  a “ 0.8X.  Again  denotes  the  pattern  due  to  Hg  in  the  y = 0 
plane,  and  G denotes  the  pattern  due  to  H in  the  x “ 0 plane. 

XX  X 

Figure  5 shows  the  equivalent  magnetic  current  for  an  open-ended 
square  waveguide  of  side  length  a radiating  into  half  space,  for  two  dif- 
ferent wavelengths.  Figure  5(a)  is  for  a “ 0.6X,  and  Figure  5(b)  is  for 
a “ 0.8X.  The  top  figures  show  the  magnitude  and  phase  of  at  y « 5a/12,  and 
the  bottom  figures  show  them  for  M at  x “ a/2.  Again  M is  zero  at  x *■  0 and 

X X 

X “ a,  and  is  large  near  both  y “ 0 and  y ■ a.  Figure  6 shows  the  principal 
plane  radiation  gain  patterns  for  the  same  waveguide-fed  apertures,  (a)  for 
a “ 0.6X,  and  (b)  for  a = 0.8X.  Again  Gg^  denote.,  the  patteim  due  to  Hg  in  the 
y “ 0 plane,  and  G denotes  the  pattern  due  to  H in  the  x ■ 0 plane. 

XX  X 

Finally,  Figures  7 and  8 show  plots  of  the  equivalent  aperture  admit- 
tance seen  by  the  dominant  mode  for  an  open-ended  waveguide  radiating  into 
half  space.  This  aperture  admittance  is  defined  by 


ap 


1 - r 

Oy 

1 + r 0 

o 


(69) 


where  F is  the  reflection  coefficient  and  Y is  the  characteristic  wave 
o o 
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impedance,  both  for  the  dominant  mode.  Our  computations  are  compared  to 
some  previously  obtained  by  Cohen,  Crowley,  and  Levis  [4].  Figure  7 shows 
the  results  for  a rectangular  waveguide  of  dimensions  a/b  = 2.25,  and 
Figure  8 shows  the  results  for  a square  waveguide  of  side  length  a.  Measured 
results  were  found  only  for  a square  waveguide,  and  these  are  also  shown  in 
Figure  8. 

IX.  .DISCUSSION 

A general  purpose  computer  program  has  been  developed  for  rectangular 
waveguides  radiating  into  half  space  through  a rectangular  aperture.  When  the 
aperture  dimensions  are  small  compared  to  the  waveguide  dimensions,  many  wave- 
guide modes  may  be  required  to  accurately  obtain  the  waveguide  admittance  matrix. 
In  this  case  it  would  be  advantageous  to  have  an  analytic  approximation  to  the 
sum  of  higher-order  mode  contributions.  However,  we  have  not  investigated  this 
possibility.  Our  program  appears  to  give  accurate  results  even  for  relatively 
small  apertures,  although  large  numbers  of  higher-order  modes  may  be  required. 

The  numerical  examples  given  in  Section  VIII  serve  to  illustrate  the 
types  of  computations  that  can  be  made  with  the  general  program.  The  program 
is  written  so  that  the  excitation  may  be  any  mode  desired,  not  necessarily  the 
dominant  mode.  The  aperture  can  be  located  anywhere  within  the  waveguide  cross 
section.  The  principal  limitation  to  the  use  of  the  program  is  one  set  by  the 
cross-sectional  size  (in  square  wavelengths)  of  the  waveguide  and  aperture. 

As  with  most  moment  soliitions,  when  the  size  of  the  aperture  becomes  too  large 
then  too  many  expansion  functions  are  required  for  a solution.  As  a rule  of 
thumb,  for  reasonable  accuracy  we  need  at  least  5 expansion  functions  per 
wavelength  for  each  component  of  current,  or  50  expansion  functions  per  square 
wavelength.  Hence,  even  on  large  computers,  one  is  limited  to  apertures  of 
the  order  of  a few  square  wavelengths  in  size, 

[4]  M.  Cohen,  T.  Crowley,  K,  Levis,  "The  Aperture  Admittance  of  a 

Rectangular  Waveguide  Radiating  into  Half-Space,"  Antenna  Lab. 

Rept.  ac  21114  S.R.  No,  22,  Ohio  State  University,  1953. 
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PART  Tim 


COMPUTER  PROGRAMS 

I . DESCRIPTION  OF  THE  MAIN  PROGRAM 

The  main  program  computes  the  complex  coefficients  V^,  which  deter- 
mine the  magnetic  current  M according  to  (1) , the  amplitudes  Pj  of  (7) 
of  the  reflected  waveguide  modes,  the  equivalent  aperture  admittance 
[1,  Eq.  (67)]  seen  by  the  Incident  mode,  the  complex  power  flowing  out 
of  the  aperture,  and  four  gain  patterns.  The  four  gain  patterns  are 
written  on  the  first  record  of  direct  access  data  set  6 so  that  they 
may  be  plotted  by  the  program  on  pages  43  and  44  of  [2].  The  main 
program  calls  the  subroutines  AY,  YMAT,  PLANE,  DECOMP  and  SOLVE  which 
are  described  in  Sections  II,  III,  and  IV, 

The  data  cards  are  read  early  in  the  main  program  according  to 

READ(1,11)  LX,  LY,  LM,  LN,  LI,  NTH,  DX,  DY,  AL,  BL,  XI,  Yl,  ER 
11  FORMAT (613,  3E14.7/4E14.7) 

The  variable*  LX,  LY,  DX,  DY,  AL,  BL,  XI,  and  Yl  are  respectively  L , 

3C 

L , &x/X,  4y/X,  a/X,  b/X,  x./X  and  y. /X  where  X is  the  free  space 
y XX 

wavelength.  See  Fig.  1.  The  variables  LM  and  LN  are  respectively  L^ 
and  L^  appearing  in  (16)  and  (17) . We  require  that 

LX  > 2 

LY  2:  1 . 

The  excitation  of  the  waveguide  is  the  LIth  waveguide  mode  where,  in 
TV 

the  program,  e^  of  (16)  is  called  the  (m+n(L  +l))th  waveguide 

mode  and  of  (17)  is  called  the  (L  +L  (L  +1)  + m + (n-l)L  )th 

m m n m m 

waveguide  mode.  ER  is  the  relative  dielectric  constant  of  (31)  and 

(32)  inside  the  waveguide.  The  four  gain  patterns  are  generated  by 

evaluating  the  plane  wave  measurement  vectors  [2,  Eqs.  (55)  tp  (58)]  at 

angles  (6  or  (j>)  equal  to  (J-1) *180 . / (NTH-1)  degrees,  J-1,2, . . .NTH. 
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PART  T\W 


COMPUTER  PROGRAMS 

I.  DESCRIPTION  OF  THE  MAIN  PROGRAM 

The  main  program  computes  the  complex  coefficients  which  deter- 
mine the  magnetic  current  M according  to  (1) , the  amplitudes  Pj  of  (7) 
of  the  reflected  waveguide  modes,  the  equivalent  aperture  admittance 
(1,  Eq.  (67)]  seen  by  the  Incident  mode,  the  complex  power  flowing  out 
of  the  aperture,  and  four  gain  pattenvs.  The  four  gain  patterns  are 
written  on  the  first  record  of  direct  access  data  set  6 so  that  they 
may  be  plotted  by  the  program  on  pages  43  and  44  of  [2].  The  main 
program  calls  the  subroutines  AY,  YMAT,  PLANE,  DECOMP  and  SOLVE  which 
are  described  in  Sections  II,  III,  and  IV, 

The  data  cards  are  read  early  in  the  main  program  according  to 

RFAD(1,11)  LX,  LY,  LM,  LN,  LI,  NTH,  DX,  DY,  AL,  BL,  XI,  Yl,  ER 
11  F0RMAT(6I3,  3E14.7/4E14.7) 

The  variables  LX,  LY,  DX,  DY,  AL,  BL,  XI,  and  Yl  are  respectively  L , 

L , Ax/X,  Ay/X,  a/X,  b/X,  x./X  and  y. /X  where  X is  the  free  space 
y J.  X 

wavelength.  See  Fig.  1.  The  variables  LM  and  LN  are  respectively  L^ 
and  L^  appearing  in  (16)  and  (17) . We  require  that 

LX  > 2 
LY  ^ 1 . 

The  excitation  of  the  waveguide  is  the  LIth  waveguide  mode  where,  in 
*TF 

the  program,  e , „ ,,v  of  (16)  is  called  the  (m+n(L  +l))th  waveguide 
tn 

mode  and  £^(^-1)1  of  (17)  is  called  the  (L  +L  (L  +1)  + m + (n-l)L  )th 
m m n m m 

waveguide  mode.  ER  is  the  relative  dielectric  constant  of  (31)  and 

(32)  inside  the  waveguide.  The  four  gain  patterns  are  generated  by 

evaluating  the  plane  wave  measurement  vectors  [2,  Eqs.  (55)  tp  (58)]  at 

angles  (6  or  ({i)  equal  to  (J-1)*180./(NTH-1)  degrees,  J«l,2, . . .NTH. 
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Minimum  allocations  are  given  by 


COMPLEX  TI(N),  GAM(NT),  P(A*N),  YHS(N*N), 
Y(NT),  V(NT),  V(N),  YWG(N*N) 

DIMENSION  GA(4*NTH),  A(NT*N),  IPS (N) 

in  the  main  program,  by 

COMPLEX  Y(NT) 

DIMENSION  A(N*NT),  Sl(LMP),  S2(LMP),  SM(LMP*LX), 
CM(LMP*LX),  CN(LNP*LY),  SN(LNP*LY) 

in  the  subroutine  AY,  by 

COMPLEX  TC(Jl),  TX(Jl),  TY(Jl),  YXX(J2),  Y(N*N) 

in  the  subroutine  YMAT,  by 

COMPLEX  P(A*N) 

in  the  subroutine  PLANE,  by 

COMPLEX  UL(N*N) 

DIMENSION  SCL(N),  IPS(N) 

in  the  subroutine  DECOMP,  and  by 

COMPLEX  UL(N*N),  B(N),  X(N) 

DIMENSION  IPS(N) 

in  the  surboutlne  SOLVE.  Here, 

N = (I.X-1)*LY  + LX*(LY-1) 

NT  = 2*LM*LN  + UN  + LN 
UP  = LM+1 
LNP  = LN+1 
J1  = <LX+1)*(LY+1) 

J2  = MAX((LX-1)*LY,  LX*(LY-1)) 

Referring  to  (22)  to  (25) , statement  A1  stores 
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In  U1  where  the  subscript  0 in  (7)  is  being  Interpreted  as  LI.  Outer 
DO  loop  19  accumulates 


-jnli  + r^3-r  'i 


- jh 


2 


in  U2.  Statement  45  puts  hP^g  where  is  given  by  (62)  in  U2. 
From  [1,  Eqs.  (45)  and  (48)], 

II  = 1 

guide 


where  is  the  transverse  electric  field  of  the  forward  traveling 

(incident)  mode  at  z = 0.  Hence  expression  (68)  reduces  to  . 

/ab 

The  V which  is  printed  upon  exit  from  DO  loop  34  is  State- 

ment 46  puts  of  [1,  Eq,  (67)]  in  Ul. 

Concerning  (59),  Nested  DO  loops  23  and  24  accumulate 

^yhs*^* 

in  U2,  Statement  47  puts  in  Ul. 

Do  loop  26  stores  G of  (8)  in  G(K).  Statement  27  uses 
TH  " (J-1*it/(NTH-1)  radians  to  store  in  P(i  + (K-1)*N).  P®  is 

given  by  [2,  Eq.  (54  + K)]  where  9 or  4>  is  TH.  For  the  2A^Ay  ^ stored 

in  P(1  + (K-l)*N)  through  P(K*N),  DO  loop  29  accumulates  2^^  in  Ul. 
Next,  G of  (8)  is  stored  in  both  G(K)  and  GA(J  + (K-1)*NTH). 

Statement  32  writes  GA  on  the  first  record  of  data  set  6 for 
possible  input  to  the  plot  program  listed  in  [2,  pages  43-44]. 
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C LISTING  GF  THE  MAIN  PROGRAM  AND  SAMPLE  DATA 
C 

//  = XEC  WATFI  V 

//G0.FT06F001  OD  OSNAME=EEOO 3<f .REV  l,D  I SP=OLD,  UNI  T=  3330* 

//  DCP=(RECFM=VS,BLKSI ZE=2596 ,LRECL=2592 ) 

//GO.SYSIN  OD  * 

4J03  MAUTZ.TIME  = ltPAGES*40 

C 

C.  MAIN  PROGRAM 

C THIS  PROGRAM  CALLS  THE  SUBROUTINES  AY, YMAT, PLANE »DECOMP t SOLVE 
COMPLEX  TI  I 50)  ,GAM(40C)  ,P(200»  ,CONJG 
COMPLEX  YHS  (2500  ) ,Y  ( 400  J ,U2,  V { 400) , U 1,  YHG(  2500) 

DIMENSION  G(4)  ,G  All  1 68)  , A (2  500)  , I PS  (50) 

REW  INO  6 

RFAD(  1,  U)  LX,LY,LM,LN,Ll ,NTH ,DX t DY , AL, BL,Xl ,Yl,EP 

11  FORMAT  (613, 3F14.7/4E14. 7) 

WRI TE( 3,12)  LX,LY,LM,LN,LI ,NTH,OX,DY, AL,BL,Xl ,Yl , ER 

12  FORMATC  LX  LY  LM  LN  LI  NTH  ' , 5X,  *DX‘ , 12X,  *D  Y* , 12X,  • AL  • /I  X ,61  3 , 
13E  14.  7/7X,*BL'  ,12X,»  Xl«  ,12X,»  Yl»  ,12  X,'  ER*  /IX,4E14.7) 

PI=3. 141593 
P?=2.*PI 
DX=0X*P2 
DY=DY*P2 
Xl  = Xl>KP2 
Y1=Y1*P2 
AL=.5/AL 
BL*.  5/BL 

41  CALL  AY  (LX,LY,LM,LN,OX,OY,AL,BL,Xl,Yl,ER,A,Y) 

42  CALL  YMAT(LX,LY,DX,OY,YHS) 

N=(LX-l  )*LY+LX*(LY-l ) 

NTa2*LM*LN4-  LN+IM 
U2=2.*Y{LI) 

J2*l 

J5  = N*(L  I-l  ) 

DO  13  J=1,N 
J1  = J 

DO  14  1=1, NT 
V (I  )=Y  < I)*A(Jl) 

JUJUN 

14  CONT  INUF 
J3=J2 

DO  15  I=J,N 
U1=YHS(  J2) 

J4=I 

DO  16  K=l,NT 
U1=U1*-A(J4)*V(K) 

J4=J4+N 
16  CONT  INUE 
YWG(  J2)  =U1 
YWG(J3 )=U1 
J2=J2+1 
J3=J3+N 

15  CONTINUE 
J2=J2+J 
J5=J5>1 

T I(J)=U2*A(J5 ) 

13  CONTINUE 

43  CALL  DECOMP  (N,  IPS, YWG) 

44  CALL  SOLVFI N,I PS,YWG ,TI ,V) 

J1  = 0 
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U2=0. 

00  19  I=1»NT 

U 1=0. 

OQ  20  J = l ,N 
J l=J  U I 

ui=uuA(  jn*v  (j) 

20  CONTIMUE 

U2=U24Ul*CnNJG(Ul  (U 
GAMd  ) = Ul 
19  CCNTINUe 

9 5 U2=(0.  f l.»»(CONJG(U2)-2.*GAM(Ln*CONJG(VlLn  n 
GAK(Ln=GAM(Ll)-l. 

CV=PI*S0RT(P2/(AL*BL*0X*0Y1  ) 

DO  39  1=1, N 
V(  n=CV*V(  I ) 

39  CONTINUE 

WRITEJ  3,35)(V(  I»,I=l,N» 

35  FORMATS  COFF  F 1 Cl  ENT  S V DIVIDED  BY  RMS  INCIDENT  6 OVER', 

V GUIDE  CROSS  SEC  TION' /(  IX,6E  11.  9)) 

WRITF(3,33)  {GAM(I),I  = l,NT) 

33  FORMAT!  'OAMPLITUDES  OF  REFLECTED  WAVEGUIDE  MODES' / (1  X ,6  El  I .9  )J 

96  U1  = (1./376.730)*(0.«1.I*Y(L1  l + U.-GAMILU  )/(  l.+GAM!LI  )) 

WRITF(3,36)  U1 

36  FORMATCOEOUIVALENT  APERTURE  ADMITTANCE  3F  INCIDENT  MODE=  ' , 2E  11.  9) 
WRITF!3,90)  U2 

^0  FORMAT  {'0(C0MPLEX  WAVEGUIDE  POW  FR  )*ET  A=  • , 2E 1 1.9) 

U2=0. 

DO  23  1=1, N 

U1=0. 

J1=I 

on  ?9  J=1,N 
U1=U1+YHS(J1)#VU) 

J l=J  UN 
29  CONTINUF 

U2=U2fV( I )^CQNJG( Ul) 

23  CONTINUF 

97  U1  = -1./(CV*CV)^(0.,1.)*U2 
WRITF(3,39)  Ul 

39  FORMAT  ( ’OICOMPLEX  HALF  SPACE  P3WERH<ETA=' ,2Ell.9  j 
CG=DX*DY/AIMAG!U2) 

OTH=PI/(NTH-l) 

P0=iaO.  /PI 
WRITE(3,25) 

25  FORMATCO  ANGLE' ,5X, ' Gl»  ,9X,  » G2  ' ,9X  , 'G3 9X,  »G9» ) 
on  26  J = 1,NTH 
TH  = ( J-1  )*DTH 

27  CALL  PLANE!  TH, LX, LY,DX,DY,P) 

TH=TH*P8 
J1=0 
J2  = J 

DO  28  K=l,9 
U1=0. 

on  29  1 = 1, N 
J1  = JU1 

U1=UUP(  JU^V!  I ) 

2 9 CONTINUF 

H=UI>XC'JNJG(  Ul) 

G (K) =CG*H 
GA(J2)  = G(K) 

J2=J?+NTH 
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28  CONTINUE 

WRITE  I 3. 3 0)  TH,(G(1)  ,I  *l  ,<V) 

30  FORMAT  { IX,F7.2,^E11.4» 

2 6 CONTINUE 
KA=J2-NTH 

32  WRI  TF<6)  (GAU)  ,J»l,KA) 

STOP 

END 

iOAT  A 

5 I 9 0 I !<}  0. 5000000E-01  0. 5000000E-01  0 .2500000E+00 

0.5000000F-01  O.OOOOOOOE+00  O.OOOOOOOE+00  0* lOOOOOOEf 01 
$ST3P 
/♦ 

// 

PRINTED  OUTPUT 

LX  LY  LM  LN  LI  NTH  OX  DY  AL 

5 19  0 I 19  0.  50000006-01  0. 5000000E-01  0 .2500000E-5-00 

8L  XI  Yl  ER 

0.5000000F-01  0.  OOOOOOOEfOO  0.  OOOOOOOF«-O0  0. 10000006+00 
COEFFICIENTS  V DIVIDED  RY  RMS  INCIDENT  E OVER  GUIDE  CROSS  SECTION 
0.1336E  + 01-0.2699E-01  0.  2 0796+01 -0.3602  E-Ol  0. 2079 6+0 1-0. 3602F-01 
0.13366+01-0.26996-01 

AMPLITUDES  OF  REFLECTED  WAVEGUIDE  MODES 
0.5  12  IF  + 00-0.27A3E-01-0.  41726-06  0. 1 8636-07  0.2043  6-01-0.1874  6-02 
0.7  153F-06-0  .121 16-07-0.65846-06  0, 1307E-07  0.  4098E-06-0.  90806-08 
-0.3752F-02  0.34436-03-0,40986-07  0.23286-09-0.1867  6-01  0.33876-03 

EQUIVALENT  APFRTURF  ADMITTANCE  OF  INCIDENT  MODE=  0 .1 1036-03-0 .1481F-02 

(COMPLEX  WAVEGUIDF  PnW£R)*6TAs  0.9504E-01  0.12686+01 

(COMPLEX  HALF  SPACE  PnWER)+ETA=  0.9504E-01  0.12686+01 


ANGLE 
3.00 
10,00 
23  .00 

30.00 
43  .00 

50.00 

60.00 

70.00 

80.00 
90.00 

100.00 

110.00 

120.00 

130.00 

140.00 

150.00 

160.00 
17(3.00 
180.00 


G1 

0.00006+00 
0.  41446-01 
0.16256+00 
0,  3528F  + 00 
0.5945E+00 
0,  862  IF+OC 
0.11236+01 
0.  1344F+01 
0.1491F+01 
0.  1543P  + 01 
0.1491F+01 
0.  1344F  + 01 
0.11236+01 
0.  8621F+00 
0.59456+00 
0.  3528F+00 
0.16256+00 
0.  4144F-01 
0.53986-12 


G2 

0 .00006+00 
0.  OOOOF  + 00 
0.00005+00 
0.  OOOOF  + 00 
O.OOOOE+00 

c.  oonoE+00 

O.Ov  •'  ‘6+00 
O.OOUCJ  + 00 
0.00006+00 
0.  COOOF^OO 
0.00006+00 
0.  00006  + 00 
0.00006+00 
0,00006  + 00 
O.OOOOE+00 
0.  00006  + 00 
0.00006+00 
0.  OOOOE  + 00 
0.00006+00 


G3 

0.00006+00 
O.OOOOE+00 
O.OOOOE+00 
0.00006+00 
O.OOOOE+00 
0.00006+00 
0.00006+00 
0.00006+00 
O.OOOOE+00 
0.00006+00 
0.00006  + 00 
0.00006+00 
0.00006+00 
0.00006+00 
0.00006+00 
O.OOOOE+00 
O.OOOOE  + OO 
O.OOOOE+00 
O.OOOOF+OO 


G4 

0. 15306+01 
0.15316+01 
0.  15326+01 
0.1533 E+01 
0.  1536E+01 
0.15386+01 
0.  15406+01 
0.1541E+01 
0.  15436+01 
0.15436+01 
0.  1543E+01 
0.15416+01 
0.  1540F+01 
0.15386+01 
0.  15366+01 
0.1533E+01 
0. 15326+01 
0.15316+01 
0.  15306+01 
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II,  DESCRIPTION  OF  THE  SUBROUTINE  AY 

The  subroutine  AY(LX,  LY,  IK,  LN,  DX,  DY,  AL,  BL,  XI,  Yl,  ER,  A,Y) 
stores  the  submatrices  defined  by  (22)  to  (25)  In  A and  the  admittances 
defined  by  (31)  and  (32)  In  Y.  More  precisely, 

A^™  Is  stored  In  A(1  + (j-l)*N), 

A^^^  Is  etored  In  A(NX  + 1 + (J-1)*N), 

A^™  Is  stored  In  A(N*NTE  + 1 + (j-l)*N),  and 
Is  stored  In  A(NX  + N*NTE  + 1 + (j-l)*N) 

where 

NX  = (LX  - 1)*LY 
N - NX  + LX*(LY~1) 

NTE  = IK*LN  + IK  + LN. 

Also, 

TF 

-jnYj  Is  stored  In  Y(j)  and 
Ti-l 

-1tiY“  Is  stored  In  Y(NTE  + j). 

The  values  of  the  variables  In  (22)  to  (32)  are  specified  by  the  first 
eleven  arguments  of  AY, 


I 
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I 

I ! 


ARGUMENT  OF  AY 

VARIABLE  IN  (22)  TO  (32) 

LX 

L 

X 

LY 

L 

y 

LM 

L 

m 

LN 

L 

n 

DX 

kAx 

DY 

kAy 

AL 

it/  (ka) 

BL 

it/ (kb) 

XI 

kx. 

Y1 

kyi 

ER 

E 

r 

Minimum  allocations  are  given  by 
COMPLEX  Y(NT) 

DIMENSION  A(N*KT),  Sl(I^),  S2(LMP),  SM(LMP*LX) 
CM(LMP*LX),  CN(LNP*LY),  SN(LNP*LY) 


where 


N = (LX-1)*LY  + lx  * (LY-1) 
I>tP  = IW  + 1 
LNP  = LN  + 1 
NT  = 2*LM*LN  + LM  + LN  . 


Statement  32  stores  loop  10  stores 

, miTox 

ipiSZ  (— --  2°-)^  M snd 

V ab  mirAx  kc 
2a 


4TiAxAye„  sin 


m 


ab 


(- 


mifAx 

2a 


mirAx 

2a 


■) 


In  S2(m  + 1) 
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whereas  nested  DO  loops  10  and  12  store 
imi(x-  + pAx) 


In  SM(mH  + (L^  + l)*p),  and 


m7T(xj^  +(p  - |)Ax) 


In  CM  (m  + 1 + (L  + l)*p). 

ra 


In  DO  loops  10  and  12,  vTM  corresponds  to  m+1  and  JP  to  p.  Nested  DO  loops 
13  and  14  store 


mr(y  + qiy) 

sin  ^ In  SN(n  + 1 (L^  + l)*q),  and 

nrrCy  + (q  - y)Ay) 

cos  r in  CN(n  -h  1 + (L  + l)*q)  . 

b n 

In  DO  loops  13  and  14,  JN  corresponds  to  n+1  and  JQ  to  q. 

DO  loop  16  stores  (22)  to  (25)  In  A.  The  DO  loop  indices  .TN,  JM, 
JP,  and  JQ  are  equal  to  n+1,  mfl,  p,  and  q respectively.  Just  before  DO 
loop  20  is  entered. 


V 2 


. nTrAy 

rnrAy 

2b 


kb  I nirAy 
^ 2b 


.Statement  27  or  29  stores  of  (31)  in  Y.  If  the  calculated  value  of 

k , 2 

(-jp)  - is  zero,  then  statement  28  replaces 


2 "6 

(-j^)  - e by  10  e_.  Just  before  DO  loop  21  is  entered, 

(C  IT 
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— j mitAx  „ . - j-  V 

~ a\  ffii*  ) V ^ 

J ^ 2a  2b 


AirAxAye 


jn  /mi 


S” 


2a 

n;tAy  ~~  J EH^  , 


^ / 2it  .xTE  £ /22S  In  A.  Nested  DO 

Nested  DO  loops  21  and  22  store  •y^jjAy 

, ^ OH  of  (23)  in  A.  Statement  33  stores 

loops  23  and  24  store  A or  v ; 

-jnY™  o£  (32)  in  '!■  "O  1°“P  “ 4^ 

loop  26  stores  *' 
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C LISTING  OF  THE  SUBRIUTINE  AY 

r 

SIBROUTINE  AY(LX,LY,LM,LN»OX,DY,  AL,RLtXl»Yl»  ERt  A,Y) 

COMPLEX  IJ,Y{400) 

DIMENSION  A (25001  ,S1 (40» ,S2(40) »SM(AOO ) ♦ CM (^00 ) t CN ( 100 ) t SN ( 100) 
U=(0.,  1.) 

32  C=S0RT(  2.5^647«J*DX*0Y*AL*BL) 

DX5=.5*DX 
LMP=LM>1 
SK  i)  = n. 

S2(l)=.  7071068»C 
00  10  JM=1,LMP 
AM=( JM-1)*AL 
DXM=DX5*AM 

IF(JM.EQ.l)  GO  TO  U 
C1=S  INI  OXM)/DXM 
S2<  JM)  =C*C1 
S 1(  JM)  = S?IJM)*C  1*AM 

11  J1=JM 
DO  12  JP=1,LX 
C 1 = (XUJP*0X)*AH 
C2=Cl-nXM 
SMI Jl) =SIN(C1 ) 

CMIJ1)  = C0S(C?) 

J1  = JULMP 

12  CONTINUE 

10  CONTINUE 
0Y5=,5+DY 
LNP=LN+1 
00  13  JN=lfLNP 
BN=IJN-1)*BL 
0YM=0Y5*BN 
J1  = JN 

DO  14  J0=1,LY 
C1=(Y1+J0*0Y)*BN 
C2=C1-0YM 
SNIJl)  =SIN(C1  ) 

CN(Jl)  = CnsiC?) 

J1  = JULNP 
14  CONTINUE 

11  continue 

LXM=LX-l 
LYM=LY-1 
NX=LXM*LY 
MY=LX'*«LYM 
N=NX+NY 
JTE  = 0 

JTM=N+(  LNP*LMP-1) 

KTE  = 0 

KTM=LNP*LMP-1 
Cl  = . 7071068 
C2=0. 

DO  16  JN=1 ,LNP 
BN=IJN-  D^BL 
BN2=PN4 BN 

IF( JM.FO. 1)  GO  TO  19 
X=RN»'DY5 
C1=S  IN(X)/X 
C2=Cl*ri'«BN 
19  DO  20  JM=1,IMP 
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IF{  ( JNfJM)  .E0.2J  GO  TO  20 
AM=(JM-l)*AL 
C 7=BN2*-AM^AM 
C8=S0RT  (f7  ) 

C7=C  7-ER 
KTE=KTF+l 
IF(C7)  27,28,29 

27  Y (KTE)  = -U*SQRT(-C7  J 
GO  TO  30 

28  C7=l.E-6*ER 

29  Y(KTE)  =-SQRT(C7) 

30  C3=Sl  (JM)*C1/C8 
C A=S2(  JM)*C  2/C8 
J1  = JN 

00  21  JO=l,LY 
C5=C3*CN(Jl I 
J1=JH-LNP 
J2=JM 

00  22  JP=l,LXM 

JTE=JTE+1 

A(  JTE»=C5*SM(  J?) 

J2  = J2+LMP 
22  CONTINUE 
21  CONTINUE 

GO  TO  ( 31)  , LY 
Jl  = JN 

DO  23  J0=1,LYM 
C5=CA*SN(Jl  ) 

Jl=Jl+LNP 
J2  = JM 

DO  2A  JP  = l,l.X 
JTE=J1  E+1 
A(  JTF)  =C5*CM(  J2) 

J2  = J2+LMP 
2A  CONTINUE 
2 3 CONTINUE 

31  IF{  ( JN-1)*(  JM-l),EQ.O)  GO  TO  20 
KTM=KTM+1 

33  Y(KTM)  =-FR/Y(  KTE) 

C5=BN/ AM 
J l=JTF-N 
00  25  J=1,NX 
JTM=JTMf I 
J1  = JU1 

A(JTM)=-C  5*A(  Jll 
2 5 CONT  INUE 

GO  TO  ( 20)  , LY 

C5=AM/BN 

DO  ?6  J=l,NV 

JTM=JTM+1 

Jl=Jl+l 

A(JTM)=C5*A{J1) 

26  CONTINUE 
20  CONT  INUF 
16  CONTINUE 
RETURN 
FND 
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THfi  SUBROUTINES  YMAT  AND  PLANE 


I II. 


THE  SUBROUTINES  YMAT  AND  PLANE  ARE  DESCRIBED  ON 
PAGES  30  TO  A1  OF  REFERENCE  2. 


C 

r 


LISTINGS  OF  THE  SUBROUTINES  YMAT  AND  PLANE 


SUBROUTINE  YMATI LX,LY,OX,DY, Y) 

COMPLEX  U,Ul»U2,U3»UA.EX,TCI  100), TX(  100), TY{  100)  ,YXX(  100) 
DX2=l.  /(OX*OX) 

DY2  = 1./  (DY*DY  ) 

DXDY=DX*DY 
NX= (LX-1  )*LY 
NY=l  LY-l)*LX 
N=NX+NY 
LXP=LX+  1 
LYP=LY  + 1 
LXM  = LX-  I 
LYM=LY-l 
U=(  0.  , 1.  ) 

UA=.1666667*U 
JST=LX+  I 
on  15  JT=1,LY 
JST=JST+ I 
YL=(JT-U5)*DY 
YU=YL^DY 
YL2=YL»YL 
YU2=YU*YU 
Yl=  (JT-l )*DY 
Y2=Y1*YI 
no  16  JS=l,LX 
XL  = ( JS-l.  5)  *DX 


,7(2  500) 


XU=XL+OX 


XL2=XL*XL 
XU2=XU*XU 
Xl=(  JS-U*OX 
X2=XI*X1 
R2=X2+-Y2 
Rl  = SORT (R2  ) 

RUl=l.-.  5*R2 

U1=RU1+Rl*( 1.-.1666667*R2)*U 
U2=R  l-RUl*U 
U3=-.5-.5*Pl*U 
EX=COS|Pl)-U*SIN(Rl) 

JST  = JST+1 
R 5=  XL  2+  YL  2 
R6=XU2+YL2 
P 7=  XL  2+ YU2 
R8=XU2+YU2 
R 1=S0RT(R5) 

R2=S0RT  (R6) 

R3=S0RT(R7) 

R4=S0RT  (R8) 

AYL  = YL*AL0G(  ( XU<-R2)  /(  XL-*-Rl)  ) 
AYU=YU*AL0G(  (XU+R4)/  (XL+R3)  ) 
AXL=XL*ALOG(  ( YU+-R3)  /{  YL+Rl)  ) 
AXU=XU*ALP.G(  (YU+R4)/  (YL+R2)  ) 
S l=AXU- AXL  + AYU-AYL 
AYL=YL*AYL 
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AYU»YU*AYU 

AXL*XL*AXL 

AXU=XU*AXU 

S3  = Xir^>AXU-XL*AXL+YU*  AYU~YL*AYL 

XY1=XL*YL 

XY2  = XU*YL 

XY3=XL>»'YU 

XY4  = XU*YU 

S5=c3333333*(  XYA*R  4- XY  3*R  3- XY2*R  XYl*P  I » f . 1 666667*  S3 
TC(  JST)  ={  SI  *Ul  *0X0Y*U2+S5*U3<-. 3333333*  (XY4*R8-XY3*R7-XY2*R  6+ XY1*R  5 
1 )*U4)*EX 
YP1=YL*R1 
YR2=YL*R2 
YR3=YU*R3 
YR4=YU*R4 

55  = . 8333333P-l*(YR4*R8-YR3*R7-YR2*R6fYRl*R5)  + . 125*{XU2*(YR  4-YR2)-X 
IL  2*(YR3-YRl)  + XU2*AXU-XL  2*AXL) 

56  = .25*0Y*?  XU2*XU2-XL2*XL2)  + .3333333*X1*DX*(YU2*YU-YL2*YL  ) 

TX(  JST  ) = .5*(  YR4-YR3-YP  2*  YR  If  A XU- A XL)  * Ul+ X1*D  XDY*  U2f  S5*U3+ S6*U4 

TX< JST) =TX{ JST)*FX/DX 

XR1=XL*R  1 

XR2=XU* R2 

XR3=XL*R3 

XR4=XU* R4 

S5=,833  3333F-  1*(  XR4*R8-XR3*R  7-XR2*R6f  XR  1*P  5)  f . 125*  ( YU2*  (XR4-XR3)  -Y 
IL2*  (XR2-XR1  )f  YU2*AYU-YL2*AYU 

S6=.25*DX*(YU2*YU2-YL2*YL2)f  .3333333* Yl*DY*(  XU2*XU-XL2*XL) 

TY(  JST)  =.5*(XR4-XR3-XR2fXRlf  AYU-AYL)*UlfYl*UXDY*U2fS5*U3fS6*U4 
TY(JST)=TY(JST)*FX/DY 

16  CONTINUF 
15  CONTINUF 

IF(LYM)  A4,44,45 

44  Jl=LXPfl 
J2=Jlf2 
TC(Jl)  = TC(J2) 

TX(  Jl)  =-TX(  J2) 

TY(Jl)  = TY(J2) 

GO  TO  46 

45  Jl=2*LXPfl 
DO  17  JS=2,LXP 
Jl^Jlfl 
TCI  JS)  =TC(  Jl) 

TXIJS)  = TX(J  1) 

TY(  JS)  =-TY(  Jl) 

17  CONTINUE 
Jl  = l 

DO  18  JT=1,LYP 
J2=Jlf2 
TC(J1)  = TC(J2) 

TXI  Jl)  =-TX(  J2) 

TY(J1)=TY(J2) 

Jl=Jlf  LXP 

18  CONTINUF 

46  J4=LXf2 
JY=0 

DO  19  JT=2fLYP 
DO  20  JS=2,LX 
J3=J4 
J4=J4fl 
J5=J4f  I 
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JY=JY+1 

YXX(JY)  = .5*(TC(  J^)  + ( JS-.5H«TC{  J5)-(  JS-3.5)<‘TC(J3»-TX(  J5)fTX(J3)  )+D 
IX2«  (TC(  J5)-2.*TC(J4)  +rC(J3)) 

20  CONTINUE 
J4=J4+2 

19  CONTINUE 
JY=0 

DO  2A  JT»l,LY 

DO  23  JS  = ltLXM 

DO  22  JO=lfLY 

JTO  = LX M*l  A8S  (JT-JQH-1 

on  21  JP=1,LXM 

Jl=JTQfIABS (JS-JP) 

JY=JY>1 
YUYI  =YXXIJU 

21  CONTINUE 

22  CONTINUE 
JY=J  Y+NY 

2 3 CONTINUE 
2A  CONTINUE 

TF(LYH.EQ.O)  RETURN 

JA=2*LXP  + 1 

JY=0 

DO  25  JT=3,LYP 
DO  26  JS=2,LX 
JY=JY+1 
J4=J4+1 
J 3=  J 4~L  XP 

YXX( JYI =(-TC( J4) fTCI J3 )+TCIJ4+l )-TC( J3+1) )/DXDY 

26  CONTINUE 
J4=J4+2 

2 5 CONTINUE 

JY  = NX 

DO  30  JT=1,LY 

DO  29  JS=1,LXM 

no  28  J0=1,LYM 

JT0=2*{  JT-JO-l 

J2=L  XM*(  lAB  SIJTOJ-  1)  /2 

DC  2 7 JP=1,LX 

JSP  = 2*IJS-JP)  +1 

Jl=J2f(lABS (JSP)+l)/2 

J Y=J Y+1 

Y( JY)=YXX(J1) 

IF(JTP*JSP.LT.OI  Y(  JY)  =-Y(  JY) 

27  CONTINUE 

28  CONTINUE 
JY=JY+NX 

29  CONTINUE 

3 0 CONTINUF 

JY=0 

J4  = LXP+2 

DO  31  JT=2,LY 

DO  3?  J9=3,LXP 

J3=J4 

.I4=J4^1 

J5=J4+l.XP 

JY  = JY+l 

YXX(  JY)  = (“TC(  J4)fTC(J3)*TC(  J5J-TC(J5“U  ) /DXOY 
32  CONTINUE 
J4=J4^2 


43 


31  CONTINUE 
JY=N*NX 

DO  36  JT=l  ,LYM 

DO  35  JS=l,LX 

DO  34  JO=l tUY 

JTQ=2*(JT-JO)-H 

J2  = LXM*(IABS(  JTQI-l)  /2 

DO  33  JP=ltLXM 

JY=JYf  I 

JSP  = 2’<«US-JP)~l 
Jl=J2+(  1A8S(JSP)  + U /2 
Y (JY  ) = YXX(J  1 ) 

IF(  JTQ*  JSP.  LT.  0)  Y(JY)«-Y(JY) 

3 3 CONT  INUF 

34  CONTINUE 
JY=JY+NY 

35  CONTINUE 

36  CONT  INUE 
JY=0 
J4=LX+2 

DO  37  JT=2»LY 
DO  3fl  JS=2,LXP 
JY=JY4-l 
J4=J4+l 
J5=J4<-LXP 
J3  = J4*'LXP 

YXX(  JY)  =.  5*(TC(J4)*-(  jT-.5»*TC(.J5)-(JT-3,5  )*TCtJ3  )~TY  ( J5  )<-TY  ( J 3)  ) ^-0 
lY2*(TC(J5)-2.nC(J4)+TC(J3» ) 

38  CONTINUE 
J4=J4>1 

37  CONTINUE 
JY= (N+1  )*NX 

on  42  JT=1,LYM 
DO  41  JS=l,LX 
no  40  JO=ItLYM 
JTQ=LX*IABS  UT-JO)+l 
DO  39  JP=l,LX 
J1  = JTQ+IABS  (JS-JP  ) 

JY=JY+l 
Y(JY)=YXX{J1  ) 

39  CONTINUE 

40  CONTINUE 
JY=JY+NX 

4 1 CONTINUE 
42  CONTINUE 

RETURN 

END 

SUBROUT  INE  PL ANE( TH» LX f LY,DX ,DY ,P ) 

COMPLEX  UtUlfP(200) 

U = (0.,  I .) 

LXM  = LX-  I 
LYM=LV-1 
NX=LXM*LY 
N=NX  4-LYM^LX 
N4=N*4 
DO  89  J=l,N4 
P( J>  =0. 

89  CONTINUE 
SN=SIN(  TH) 

CS=COS  (TH» 
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X2=0X*CS 
X3=.5*X2 
Sl=-SIN(X3)  /X3 
S2=SI+S1*SN 
DO  81  JP«l»LXM 
S5  = JP*X  2 

Ul=S2*(COS(  S5)  + U*SIN(S5)) 
Jl=JP 

DO  87  JO=ltLY 
P{J1)=U1 
J l=J  1 + L XM 
87  CONTINUF 

81  CONTINUE 
IFILYM.fiQ.O  ) GO  TO  90 
DO  82  JP*1»LX 

S5  = ( JP-.5)’*'X2 
Ul=Sl*{COS(  S5)  + U*SIN(S5)I 
J1=N+NX+JP 
DO  88  JQ=l,LYM 
P(J1 )=U1 
J 1=J  l+L  X 
8 8 CONTINUF 

82  CONTINUE 

90  Y2=DY*CS 
Y3=.  5*Y2 
S1=-SIN(Y3 I/Y3 
S2=Sl*Sl^SN 
Jl=2*N+NX 

IF(LYM,FO.O)  GO  TO  91 
DC  8'.  J0=1,LYM 
S5=JO*Y2 

UI=S2*1CCS(S5  »+U*SIN(S5H 
DO  8A  JP=lfLX 
Jl=Jl+l 
P(J  l)  = Ul 

84  CONTINUF 

83  CONTINUF 

91  DC  85  JO=ltLY 
S5=IJ0-.5)*Y2 
Ul=Sl*(CCS  (S5  )+U*SYN(S5n 
DO  86  JP=ltLXM 

Jl=Jl+l 
P(J  n = ui 
86  CONTINUF 

85  CONTINUF 
RETURN 
FND 


45 


IV,  DESCRIPTION  OF  THE  SUBROUTINES  DECOMP  AND  SOLVE 


The  subroutines  DECOMP (N»  IPS,  UL)  and  SOLVE (N, IPS, UL,B,X)  use  the 

method  of  Gaussian  elimination  and  LU  decomposition  described  in  [5 , 

Section  9]  to  solve  a linear  system  of  equations  with  complex  coefficients. 

DECOMP  and  SOLVE,  based  on  the  FORTRAN  programs  on  pages  68  and  69  of  [5], 
n3 

require  roughly  multiplicative  operations  to  solve  a system  of  n linear 
equations  tdiertsas  the  subroutine  LINEQ  on  pages  28  and  29  of  [2]  requires 

3 

roughly  n multiplicative  operations  to  solve  the  same  system. 

The  input  into  DECOMP  is  N and  the  matrix  [A]  of  coefficients  of  the 
set 

[A]x  " ? (69) 

of  N linear  equations  stored  by  columns  in  UL.  N ^ 2.  The  output  from 
DECOMP  is  IPS  and  UL.  DECOMP  does  not  change  N.  SOLVE  uses  N,  IPS,  UL, 
and  the  elements  of  b stored  in  B to  calculate  and  store  in  X the  elements 
of  the  solution  vector  x.  SOLVE  does  not  change  any  of  the  input  variables 
N,  IPS,  UL  and  B.  Minimum  allocations  are  given  by 

COMPLEX  UL(N*N) 

DIMENSION  SCL(N),  IPS(N) 

in  DRCOMP  and  by 

COMPLEX  UL(N*N),  B(N) , X(N) 

DIMENSION  IPS(N) 

in  SOLVE. 

Assuming  that  the  reader  is  familiar  with  (5,  Section  9],  we  will 

attempt  to  explain  interchanging  rows  and  scaling  by  rows  [5 , Sections  10 

and  11] . Interchanging  rows  is  necessary  to  avoid  pivots  which  are  close 

(2)  (3) 

to  zero.  The  divisors  8^2  and  a^g  of  [5,  Section  9]  are  called 
pivots.  Scaling  the  matrix  by  rows  facilitates  the  search  for  pivots  by 
insuring  that  the  number  with  the  largest  magnitude  is  most  likely  the 
best  pivot. 

It  will  be  shown  that  interchanging  the  ith  and  Jth  rows  of  th^ 
array  UL  in  storage  at  any  stage  of  the  computation  is  equivalent  to 

[5]  G.  E.  Forsythe  and  C.  B.  Moler,  "Computer  Solution  of  Linear  Albegraic 
Systems,"  Prentice-Hall,  Inc,,  Englewood  Cliffs,  Hew  Jersey,  1967. 
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Interchanginp  the  ith  and  jth  elements  of  b.  The  above  statement  Is 
certainly  true  If  the  row  interchange  occurs  at  the  very  beginning  when 
[A]  resides  in  UL.  After  zeros  have  been  introduced  in  the  below  the 
main  diagonal  positidns  of  the  first  m columns  of  A, 


UL  - 


a^2- 

a^2. 


(70) 


11  11 
where  A is  an  m by  m submatrix.  The  elements  of  A on  and  above  the 

12  22 

main  diagonal  and  all  the  elements  of  A and  A are  elements  of 

[MM  , ...  M,A]  where  M.  is  defined  in  [5,  Section  9].  The  elements  of 
m—1  1 j 2^ 

A strictly  below  the  main  diagonal  and  all  the  elements  of  A are  ele- 
ments of  [M  M , , . , M, ] . Hence, 

m m-1  1 ’ 


^ B^^  o' 

■ A^2 • 

1—i 

i 

. 0 A^^ . 

P 


(71) 


where  is  a lower  triangular  matrix  whose  main  diagonal  elements  are 
equal  to  unity  and  whose  elements  below  the  main  diagonal  are  equal  to 
those  of  A^^,  C is  an  upper  triangular  matrix  whose  elements  on  and 
above  the  main  diagonal  are  equal  to  those  of  A'*'  , U is  an  identity 
matrix,  and  b'  is  b with  its  elements  permuted  in  accordance  with  any 
previous  row  interchanges  which  may  have  occurred.  From  (71), 


.21,12  . .22  , 

A A + A J 


X =•  ?’ 


(72) 


At  the  stage  ' • it-  »iiven  by  (70),  the  row  Interchange  is  limited 

21  22 

to  the  last  mfl  ro«-^  -■>?  ti.  hence  only  affects  A and  A . From  (72), 

7.1  22 

an  Interchange  of  two  rows  of  A and  A is  equivalent  to  interchanging 
the  corresponding  two  elements  of  b . 

In  the  rovr  Interchanges  are  not  dona  physically,  but  mentally 

by  means  of  the  variable  IPS.  The  Ith  row  of  the  row  Interchanged  Array 
is  the  IPS(I)th  row  stored  in  UL. 
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Scaling  by  rows  consists  of  multiplying  all  the  elements  of  the  jth 
row  of  A,  j=l,2,...N,  by  the  reciprocal,  called  SCL(j)  In  DECOMP,  of  the 
magnitude  of  the  largest  element  In  the  jth  row.  In  DECOMP,  the  rows  of  A 
are  not  scaled  physically.  The  effect  of  physical  scaling  would,  as  de- 
duced from  the  algorithm  In  [3,  Section  9],  be  to  multiply  any  jth  row 
element  involved  in  a search  for  a pivot  by  SCL(j).  DECOMP  obtains  the 
effect  of  scaling  by  Incorporating  SCL  into  the  searches  for  pivots, 

A verbal  flow  chart  of  DECOMP (N, IP S,UL)  is  as  follows.  DO  loop  5 
Initializes  IPS  and  obtains  SCL.  In  DO  loop  2,  the  operation  consisting 
of  the  sum  of  the*  abcolute  values  of  the  real  and  imaginary  parts  ought  to 
be  more  efficient  than  the  usual  magnitude  Involving  the  square  root  of 
the  sums  of  the  squares  of  the  real  and  imaginary  parts. 

DO  loop  17  performs  the  premultiplication  by  Mj^  where  M^^  is  defined 
in  [5,  Section  9].  DO  loop  11  decides  that  the  (IPS(IPV) ,K)th  element  of 
IIL  will  be  the  pivot.  If  we  were  using  the  method  of  physical  row  inter- 
changes we  would,  at  this  point.  Interchange  the  Kth  and  IPVth  rows 
physically.  Statements  14  and  the  two  statements  following  it  carry  out 
the  corresponding  mental  row  Interchange  by  interchanging  IPS(IPV)  and 
IPS(K).  Nested  DO  loops  16  subtract  m^j^  times  the  IPS(K)th  row  of  UL 
from  the  last  K+1  columns  of  the  IPS(I)th  row  of  UL.  See  [5,  Section  9] 
for  the  definition  of  Note  that  m^j^  is  such  that  the  previous  row 

operation  extended  to  the  Kth  column  would  annihilate  the  (IPS(I),K)th 
element  of  UL.  Statement  18  stores  m^j^  in  the  (IPS(l),K)th  position  of 
UL. 

The  subroutine  SOLVE(N,IPS,UL,B,X)  obtains  the  solution  x to 

[L][U]x-S:'  (73) 


where  [L]  is  a lower  triangular  matrix  whose  diagonal  elements  are  all 
unity  and  [U]  is  an  upper  triangular  matrix.  Here,  U^^ , j ^ 1,  is  stored 
in  the  (IPS(l),  j)th  position  of  the  input  array  UL.  Also,  “ 1 and 
L^j,  i > j is  stored  in  the  (I?S(l),j)th  position  of  UL.  The  ith  element 
of  b*  is  stored  lu  B(IPS(1)),  As  ia  [5,  Section  9],  we  let 


1 


(74) 


(U]x  “ 


->• 

y 


such  that  (73)  becomes 


[L]y  » (75) 

The  solution  to  (75)  la  given  by 

yi  = (76) 

i-1 

y.  = - I L y,  , i-2,3,...N 

j“l  •* 

The  solution  to  (74)  Is  given  by 


% 

'nn 


■ ).L 

’■i — — 


(77) 


1»N-1,  N-2,...l 


A verbal  flow  chart  of  the  subroutine  S0LVE(N,1PS,UL,B,X)  Is  as 
follows.  DO  loop  2 stores  of  (76)  In  X(I).  The  Index  J of  inner 
DO  loop  1 Is  the  surama:  ion  index  j appearing  In  (76) , DO  loop  4 stores 
x^  of  (77)  in  X(i)  for  1 * ”+l  - IBACK.  The  Index  J of  inner  DO  loop  3 
is  the  summation  index  j appearing  in  (77). 
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C LISTINGS  OF  THE  SUBROUTINES  OFCOMP  AND  SOLVE 
C 

SUBRQJT  IMC  OFCOMPIN,  IPStUL » 

COMPL<=X  UL(  ?500)  ,PI  VOT,EM 
niMFMS  In^'  SCI  (50),  IPS(50) 

•)P  5 1 =1,N 
I PS ( I)  = I 
RN  = 0. 

Jl=  I 

DO  2 J=1,N 

JLM=ABS  (PEAL  (IJL  (Jin  1 fABS  ( A IMAG(  UL  ( J U)  ) 

Jl  = Jl  + N 
IF(RN-JL^) 

1 RN=ULM 
7 CONTINUE 
SCI  ( I)  = 1.  /RN 
5 CONTINUE 

NML=m-  1 
K2  = n 

on  17  K=l,MMl 
RIG=n. 
on  11  I =K,N 
IP=  IPS(  I ) 

I PK  = I P+  K? 

5IZE=(  ABS(RSAL(UL(  IPK)  ) )fABS(AI  MAG(UL(1  PK)  n )*SCl  (IP) 
IF(SI7E-BIG)  11,11,10 
10  R1G=SI2E 
I PV  = I 

1 1 CONT  INUE 

IF(IPV-K)  LA, 15, 14 

14  J=  IPS(K  ) 

IPS(K)=IPS(IPV) 

IPS( IPV  ) = J 

15  KPP=t PS (K) +K? 

P 1V0T=JL(KPP) 

KP1=K«-1 

no  16  1=KPI,N 
KP=KPP 

IP=  IPS(  I )<-K  ? 

FM=-UL(  I P)/  PIVOT 
Ifi  UL ( IP  ) = -PM 
DO  16  J=KP1  ,N 
IP= IP+N 
KP  = KP+-N 

UL  ( IP  ) = UL(  IP)4-EM*UL(  KP) 

16  CONTINUE. 

K2=K2+N 

17  CrNTINlJC: 

PFTUPN 

END 

SUBRQJT  INF  SOLVEIN, IPS ,UL,B,X) 

COMPl  FX  UL(  2500)  ,B(50)  ,X(50)  , SUM 
DIMENSIPN  IPS  (50) 

NP  1=N4-  1 
I P=  I PS  ( 1 ) 

X(  1)  =B(  IP) 
nn  ? i=2,N 
IP  = IPS(  I) 

IPB-  IP 

I M 1=1 -I 
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SUM=0. 

DO  I J=MMl 
SUM=SUM4-UL(  IP)*X  (J  ) 

1 IP-IP+N 

2 X(i)=B{  1PB)-SIJM 
K 2=N*(N-1) 

IP=IPS(N)+K2 
X('J)  = X(N»  /UUIP) 

on  ^ iBArK=?,N 

I=V)P  l-I  BACK 

K2  = K2-N 

IPI  = IPS(I  )<-K2 

IPl  = I+l 

SUM=0. 

IP=IPI 

on  3 J = IP  ItN 
I P= I P+N 

3  SIJM=SUM+UL(  IP)*X(  J) 

^ X{  I ) = (X  { n-SUM)/UL(IPU 
RETURN 
END 


i 

I 

1 

! 

I 
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